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Abstract 

Molecular dynamics is usually performed using specialized software packages. 
We present methods to perform exploratory studies of various aspects of molecular 
dynamics using MATLAB. Such methods are not suitable for large scale applica- 
tions, but they can be used for developement and testing of new types of interactions 
and other aspects of the simulations, or simply for illustration of the principles for 
instruction and education purposes. One advantage of MATLAB is the integrated 
visualization environment. We also present in detail exploration of forces obtained 
from 3-body potentials in Molecular Dynamics in this framework. The methods are 
based on use of matrices and multidimensional arrays for which MATLAB has a 
set of both linear algebra based as well as element-wise operations. The numerical 
computations are based on these operations on the array objects rather than on the 
traditional loop-based approach to arrays . The discussed approach has a straight- 
forward implementation for pair potentials. Applications to three-body interactions 
are the main aspect of this work, but extension to any general form of n-body inter- 
actions is also discussed. The methods discussed can be also applied without any 
change using the latest versions of the package GNU OCTAVE as a replacement 
for MATLAB. The code examples are listed in some detail, a full package of the 
MATLAB and OCTAVE codes is available for download. 

1 Introduction 

Modelling of atomic motions in terms of classical Newtonian equations is usually called 
molecular dynamics (MD). Such modelling is considered useful in many areas of material 
research, chemistry, molecular biology, with various justifications for the applied approx- 
imations. In most applications large ensembles of atoms are treated using specialized 
software. In this paper we describe an alternative approach and describe methods which 
can be useful in design and exploration of new interactions, or simply in education. While 
the large systems require optimization of the computations and often possibility of paral- 
lel computation, the presented MD models are supposed to be of the size up to tenths or 
possibly hundreds of particles. We show that for such relatively small systems of model 
atoms the computations can be carried out inside the systems MATLAB or OCTAVE. 
These systems should be relatively well known, but we give detailed references and small 
presentation of features relevant for this work in the Appendix El 

It is generally well known that MATLAB and related systems have interpreted pro- 
gramming languages and this can make execution of large programs rather slow. On the 
other hand, they have the ability to perform mathematical operations on complicated 
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structured objects in one single statement using ordinary mathematical notation. This 
property is essential for our approach. A simple MD simulation code with pair potentials 
can be written in about 30 statements. These can be complemented by another set of 
routines of a comparable size (about 20 to 30 statements) for quite powerful real time 
visualizations. See Appendix |A] for some special aspects in MATLAB language. 

The system MATLAB and its open source relative GNU OCTAVE have emerged as 
"matrix laboratory", i.e. a system for working with matrices. In their latest versions they 
have been extended to handle not only matrices and vectors, but also multidimensional 
arrays. The systems are based on the use of an internal programming language. The 
coding is made easier by an internal on-line manual system, which explains in detail the 
various built in functions. The disadvantage of computing with MATLAB is that the 
language is interpreted and not compiled, thus not efficient in large scale calculations. 
On the other hand, the internal operations are based on optimized FORTRAN and C 
libraries and thus perform very well and very efficiently. Most of the statements about 
MATLAB are valid to OCTAVE, except the graphical cabilities. MATLAB has a very 
advanced graphical system, its counterpart in OCTAVE has long been slightly different 
and much simpler, but in the latest versions it is attempted to become very compatible. 

For exploratory work and small projects the fact that the programs are interpreted 
and not compiled is of a great advantage. Codes can be quickly adapted to any new 
required tasks and tested in matter of seconds, since no compilation and linking cycles 
are required. The codes discussed in this paper and additional material are available for 
download [1]. 



2 Molecular dynamics and Verlet algorithm 

Generally, molecular dynamics refers to classical mechanics modeling of atomic systems 
where the atoms are represented by mass points interacting with their surroundings, 
including the other atoms, by forces which can be represented by potentials, schematically, 
for each atom there is a Newton equation 

m,f-- = -Vy(....,r,,....) (1) 

In many cases the coupled second order equations are solved using the well known Verlet 
algorithm. This is a special difference scheme for solving Newton-type equations 

^ r{t) = a{t) (2) 

using the following difference scheme 

r(t„+i) = 2r(t„) - r{tn-i) + a{tn)At^ + O(At^) (3) 

We illustrate the use by the code for one particle moving in attractive Coulomb potential, 
with the possibility of changing it to the so called soft Coulomb, where the 1/r form is 
modified. Such code is very popular for teaching demonstrations, but in our context it 
is useful to test the accuracy and stability of the solutions in this very simple case. The 
code (we leave out the simple input and definition statements) is 

1 . t=0 : tmax/nt ime : tmax ; 

2. r=zeros (2 ,ntime) ; 
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3. r(: ,1)= rO; 

4. r(: ,2)= rl ; 

5. for nt = 2:ntime 

6. Rt=r(:,nt); accel=-Cm * Rt/sqrt(a"2+sum(Rt . "2) ) "3; 

7. r(:,nt+l)= 2*r(:,nt)- r(:,nt-l) + accel *( t (nt) -t (nt-1) )"2; 

8. end 

Here in line 6 we see that the vector Rt representing r{t) is assigned the latest known 
position and the accelerations vector is evaluated according to 

r 1 

Line 7 implements the Verlet algorithm as given above in eq. [31 This line of code with 
only occasional small modifications is used in all the discussed applications. Note that 
this operation is performed on a vector. Later, when we consider more than one particle, 
the symbol r can represent a "vector of vectors", then this line will become 

r(:,:,nt+l)= 2*r(:,:,nt)- r(:,:,nt-l) + accel *( t (nt) -t (nt-1) )~2; 

Similar notation is available in the new versions of FORTRAN. 

The simple code for this classical central force problem is very illustrative and useful 
for exploring the stability and precission of the numerical method. When the parameter a 
is zero, we have the Coulomb-Kepler problem and the motion is periodic along an ellipse. 
If the a is different from zero, the motion is generally not periodic and we get a rosette- 
like motion. This elementary example is included because the codes in this paper use the 
same Verlet implementation. 



3 Two body interactions 



The two-body forces considered here are forces between mass points, so that the potential 
energies can only be functions of the distances between the particles, perhaps with addition 
of global forces of the environment on each particle, but we shall limit this discussion to 
the pair forces only, a global force is analogue to the first example of planetary motion 
above. 

The total potential energy is given by 

Vtotiru ,r^) = J]y(|ri-r^|) (4) 

i<j 

and the force on each particle is obtained by gradient operation 



^ / \ sr^ dV(u) 

TUkYk = -VkVtot{ri, = 

The symbol 



u-R 

denotes the (m, n) element of the matrix obtained by applying function 

dV{u) 
du 



(5) 



(6) 
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on the matrix R. We will generalize this notation to any type of objects, e.g. 

dV iu) 



du 



(7) 



u=Q 



ijk 



denotes the (i, j, k) element of the corresponding 3-index array. Without the indices, the 
discussed symbol represents the whole matrix or array. 



4 MATLAB implementation 

We shall now show in detail how solutions of eq. [5] can be implemented in very short 
codes in MATLAB. 

For one- dimensional world, to generate an object holding all the distances dij = Xj — 
Xi, we can proceed as follows. We can aim at a matrix D with elements dij. In traditional 
programming such a construction would require two loops. In MATLAB this can be done 
in one statement. Provided that the positions are in a vector x = {xi,X2 xn) 

D=ones(N,l)*x-x'*ones(l,N) ; 

In mathematical notation this means 

[ Xi X2 X3 ... Xn ] 



Xi 
X2 
X3 



Xn 



1 1 1 



and the following matrix operation 





Xi 


X2 


Xs ■ ■ 


■ Xn 




Xl 


Xl 


Xl ■ ■ 


■ Xl 




Xi 


X2 


X3 ■ ■ 


■ Xn 




X2 


X2 


X2 ■ ■ 


■ X2 
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Xi 


X2 


Xs ■ ■ 


■ Xn 




xs 


xs 


X3 ■ ■ 


■ Xs 




Xi 


X2 


X3 ■ ■ 


■ Xn 




Xn 


Xn 


Xn ■ ■ 


■ Xn 



(9) 



results in the matrix of all scalar distances whose elements are the distances 



The diagonal of the matrix D is clearly zero, and that could cause some problems in 
inversions and similar operations. Besides, the forces "on itself are not considered, so 
that there should be a mechanism of keeping the diagonal out of the calculations. The 
simple mechanism in MATLAB is to use "not" operation on the unity matrix, given by 
the built-in function eye(N) 

U= eye(N) 

The matrix Un is N x N matrix with one everywhere except on diagonal, obtained by 
Un = ~eye(N) 
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In MATLAB there are two types of multiplication: the mathematical product of matrices 
A * B 

and the element by element, {A ■ *B)ij — A^Bj, appearing in code as 
A.*B 

It should be rather clear that this can be used as a decomposition by projection operators 
into a part which has zeros on the diagonal and the diagonal part itself: 

C= Un .* A + U .* A; C= (~U) .* A + U.*A; C= ~eye(N).*A + eye(N).*A 

In the above code, the three assignments perform the same operation, assigning the matrix 
A to C, but the first term in the assignment is always the non-diagonal part. This 
decomposition can be used in various combinations, as we shall see below. 

The code for our molecular dynamics engine has only one disadvantage: for speed 
it is necessary to write the function for the force explicitly, but in principle it could be 
specified by a call to function script. Here we shall use the explicit form. 

The whole code for the time development of the system of particles, i.e. the whole 
molecular dynamics engine in this approach is 

1. Umat=zeros(2,N,N) ; Umat (1 , : , : )=~eye(N,N) ; Umat(2, : , : )=~eye(N,N) ; 

2. x=zeros(l,N) ; y=x; F=zeros(2,N) ; RR=zeros(2,N,N) ; 

3. N=input('No of particles > '); ntime=input( 'No of time steps > 

4. R=zeros(2,N,ntime) ; t=0:dtime:ntime*dtime; 

5. start_solutions (R,N,iitime) ; 
6. 

7. for m = 2:ntime 

8. x(:)=R(l, : ,m) ; y( : )=R(2, : ,ni) ; 

9. dx= ones(N,l)*x - x' * ones(l,N) ; 

10. dy= ones(N,l)*y - y' * onesd.N) ; 

11. dr=sqrt(dx. "2 + dy."2) + eye(N); % eye(N) avoids div. by zero 

12. RR(1, : , :)=dr; RR(2, : , : )=dr ; 

13. dR(l, : , :)=dx; dR(2, : , : )=dy; 

14. forces = (-a*A*exp(-a*RR) +b*B*exp(-b*RR) ).* ( - dR ./ RR) ; 

15. F( : , : )=STim(f orces . *Umat ,3) ; 

16. R(:,:,m+1)= 2*R(:,:,m)- R(:,:,m-1) + F/M *( t(m)-t(m-l) )~2; 

17. end 

Note that the above is not a pseudocode, but actual code. The fifth line is a shorthand 
notation for starting the solution. The Verlet method must be started by assigning the 
starting time values of positions and the volues at the first time step. Otherwise, the first 
5 lines effectively declare the variables. The 3-index array of positions storing the whole 
simulation is R - for all times, the value of x and y (and z in 3-dimensional case) are the 
present positions, evaluated in previous step, are equivalently written as 

R(l:2,l:N,l:ntime) x(l:N) y(l:N) 

i.e. specifying the whole definition range gives the object itself. We have chosen to present 
the case of particles moving in 2-dimensional space. The 3-dimensional version of this 
code is obtained by only adding relevant analogues for z, dz and changing the value 2 
to 3 in all the zeros (2, . . . .) statements. It is naturally possible to write the code so 
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that it includes both two and three-dimensional versions, but we present here the more 
transparent explicit form used in the listing. The actual loop over the time mesh is 
shown in lines 7-17. Line 8 assigns the present time values. The matrices dx and dy are 
discussed above, they hold the components of the distances between the particles. We use 
multidimensional arrays to perform all the operations, therefore there is a new complexity 
element. Since the operations on the forces components are evaluated in one statement, 
line 14, the dr is duplicated into the appropriately dimensioned RR(1 :2, 1 :N, 1 :N) 

For the understanding of the more compact code which we shall use further on 
we reformulate the code in the above example. First we repeat the relevant part to be 
changed 

8. x(:)=R(l, : .in) ; y( : )=R(2, : ,m) ; 

9. dx= ones(N,l)*x - x' * ones(l,N) ; 

10. dy= ones(N,l)*y - y' * ones(l,N) ; 

11. dr=sqrt(dx. "2 + dy.''2) + eye(N); % eye(N) avoids div. by zero 

12. RR(1, : , :)=dr; RR(2, : , : )=dr; 

13. dR(l, : , :)=dx; dR(2, : , : )=dy; 

14. forces = (-a*A*exp(-a*RR) +b*B*exp(-b*RR) ) . * ( - dR . / RR) ; 

and below we see the more compact code. This way of coding is now applicable for both 
two-dimensional and three-dimensional physical space. First the instantaneous array of 
vectors is assigned to r (replacing x and y). Then the vector multiplications are replaced 
by constructions explained in the appendix. 

8. r(: , :)=R(: , : ,m) ; 

9. dR=r(: , : ,ones(l,N)) - permute ( r( :,: ,ones (1 ,N) ) , [1,3,2] ); 
10. 

11. RR(1 , : , : )=sqrt( sum(dR. "2, 1 )) + "Umat ; % "Umat avoids div. by zero 

12. RR=RR([1 1] , : , :) ; 
13. 

14. forces = (-a*A*exp(-a*RR) +b*B*exp(-b*RR) ) . * ( - dR . / RR) ; 

15. F( : , : )=sum(f orces . *Umat ,3) ; 

16. R(:,:,m+1)= 2*R(:,:,m)- R(:,:,m-1) + F/M *( t(m)-t(m-l) )~2; 

17. end 

The operation permute is a generalization of the transpose operation, where the permu- 
tation of indices is specified by the second variable of the operation, i.e. if A=ones(N,N), 
both of the following calls 

permute ( A, [2,1] ) -transpose (A) 
permute ( A, [1,2] )-A 

would return zero matrices. 

In the last code there is still too much calculation going on. It can be simplified by 
having the forces predefined 

8. r(: , :)=R(: , : ,m) ; 

9. dR=r(: , : ,ones(l,N)) - permute ( r( :,: ,ones (1 ,N) ) , [1,3,2] ); 

10. rr(: , :)=sqrt( sum(dR."2,l )) + eye(N); 

11. RR(1, : , :)=rr; RR=RR( [1 1] , : , : ) ; 

12. forcesd, : , :)= (-a*A*exp(-a*rr) +b*B*exp(-b*rr) ); 

13. forces=f orces ( [ 1 1],:,:); 

14. forces = forces.* ( - dR ./ RR); 



6 



15. F( : , : )=siam(f orces . *Umat ,3) ; 

16. R(:,:,m+1)= 2*R( : , : ,m)- R(:,:,m-1) + F/M *( t(m)-t(m-l) )~2; 

17. end 

The calculation might become faster for very large number of particles but the code in 
the previous listing is more understandable. For small numbers of particles the first 
formulation might even perform faster in some cases, simply because the code is shorter. 



4.1 Two body forces reformulated 

In this part we reformulate the two-body problem in order to generalize to the n-body 
case. The total energy is a sum over the pairs 



tat — / , ^ I] 

i>j 



(10) 



which can be rewritten as, ( provided that Vij — Vji which is always satisfied ) 

^tot^^EE^^^- (11) 

and the condition j ^ i is automatically reflected by the properties of the interactions. 
The force on m-th particle is obtained by 



v^yt 



m 'tot 



(12) 



In evaluating each of the contributions in the sum 

dV{u) 



du 



(13) 



we must realize that (note that this is the simplest notation, already including symmetry) 



(14) 



Due to the properties of the ^-symbols this can be written as 



dV 
dT 



(15) 



where represents the unit vector in the direction Tij and the terms after summation 
will become e^j and Cj^. Since the index m appears at different positions, after the 
summations are performed, the two terms will correspond to two different positions of the 
matrix. For simplification of the notation we will just consider the notation e = {X, Y, Z) 
for the x,y,z components of the unit vectors and consider now just the x-components Xij; 
the whole matrix will be then denoted by X. Thus 



d 



E 



dV{u) 



du 



dV iu) 



du 



(16) 
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We observe that the right hand side of the last eq. 
in the first term we use the fact that 



dV(u) 



du 



dViu) 



du 



is in the form of matrix product, 



dV 



dr 



jm 



smce r. 



mj 



and the two terms can be brought to the same form as the second, 



d ^ 



X 



dV 
dr 



(17) 



i.e. the matrix product of matrix X and the derivative matrix, not element-wise product 
encountered before. The matrix notation conveniently denotes the result, but it is not 
needed for the practical evaluation (it takes times longer time). The element-wise 
product of the two matrices (the second transposed) and summation over the second 
index will perform the same operation. In MATLAB notation. 



A= B*C; 

VEC=sum( B .* C. 



7. A matrix N x N VEC(k) = A(k,k) 



2 ); 



The reduction to a vector shown in eq. [17] can be seen as a general rule which will be 
useful below, in the three-body case. This formulation is already implemented in the 
listings in section HI 



5 Three body interactions 

Consider now a system of N mass points whose interactions are described by both pairwise 
distance dependent forces and additional three-body forces with dependence on the angles 
between the "bonds" to the neighboring points, a situation usually found in molecular 
dynamics. The two-body part of the forces can be treated as above. To evaluate quantities 
related to triplets of objects, as e.g. angles between vectors connecting the particles 
(known as bond angles), a set of all pairs would be first computed, then a loop would run 
over that set, effectively including a triple or even four-fold loop, depending on algorithms 
chosen. All further evaluations of various functions of distances and bond-angles would 
involve writing efficient loops for manipulations, evaluations and summations of energies 
and forces, and solving the differential equations for the coupled system. Here we discuss 
the extension of the approach discussed above to cover also the three-body case. We 
show the methods on the example of Stillinger and Weber three-body potentials. We 
show that also for this case the operations can be programmed in MATLAB in simple 
scalar-like notation and from the point of view of efficiency, the mathematical operations 
on matrices,or even multidimensional arrays, can be performed with speeds comparable 
to compiled and optimized code. 



5.1 Stillinger- Weber potential 

A typical three-body potential is the Stillinger and Weber potential defined (additional 
details not important for the present discussion can be found in the original reference) as 

Ksi^(ri,r2,r3,....,r7v) = ^Vadri -r^l) + ^ ^{(ri, r^, r^) (18) 

i<j i<j<k 
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and the force on each particle is still obtained by gradient operation 

mkfk = -VfcVtot(ri, , Tiv) (19) 

exactly as in eq. [5], but now the simple evaluation resulting in sum of forces from the 
individual neighbours can not be used as in the second part of eq. [5l 

The V2{rij) part contributes to the forces as described in the previous section. The 
new three body part is by Stillinger and Weber chosen as 

V3(ri, Tj, Ffc, ) = h{rij,rik, Ojik) + h{rji, rjk, Oijk) + h{rki, rtj, 9ikj) (20) 
where 6ijk is the angle between rji and r^j 

cos 9ikj = — ^ ■ T^^^^ — ^ cos 6ikj = eij ■ e^fc (21) 

\rk-ri\ \rj-ri\ 

at the position of the i-th particle. We also introduce here the notation for the unit 
vectors Sij in direction of the distance vectors r^j = rj — rj. The functions h{ri, r2, 9) are 
chosen by Stillinger and Weber as 

hijij, Tik, 9jik) = g{rij)g{rik) ^cos%fc + (22) 

where the functions g{Tij) mainly play a role of cut-off functions. We see that the most 
important feature of three-body part are the angles between the lines connecting the pairs 
of particles. This is the element which is going to be investigated in detail in the following 
section. 

For later use we denote 

9jik) (23) 

a notation useful in discussion of the evaluation of the forces below. Note that in this part 
we have followed to high degree the notation used by the authors in the original paper. 
Below we shall use our notation with stress on the pair-separability of this interaction, as 
discussed below. The main feature is the rather simple rewriting by the secongd part of 
eq. EI 



5.2 Array operations for three-body quantities 

In eq. |9]we have introduced the fast way to compute and represent the distances between 
all particles - and in the following used both the scalar distances and their vector versions. 
These quantities have two indices and are thus conveniently treated as matrices, or N x N 
objects. When we wished to keep all the space components in the same object, we ended 
with 3 X N X N array objects. For three body interactions the potential terms Vijk of eq. 
|23] are at the outset N x N x N array objects. 

To start we imagine a simplified term: instead of eq. [22] ( which we have already 
rewritten according to a more transparent notation of eq. [2T] ) 

Vijk = g{rij)g{rik) Q + 6^^- ■ 6^^^ (24) 

we first consider only the first two terms in the product 

Vijk g{rij)g{rik) (25) 

The objects with indices ijk which would represent the g{rij) and g{rik) can be built from 
the matrices D as shown schematically in figure [TJ 
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D3=zeros(N,N,N) ; 

D3(: , : , :) = D2 ( : , : ,ones(N,l)) 

This corresponds to the first cube in figure [H and it is going to be multiphed by a similar 
object obtained by permuting the indices as shown schematically in this listing: 

D3=zeros(N,N,N) ; 

D3(:,:,:) = D2 ( : , : , ones (N, 1) ) ; 
Gl=gfun(D3) ; 

G2=gfun( permute ( D3, [ 1 3 2 ] ); 
V3= Gl .* G2; 

where the two matrices are multiplied element-wise. 



Figure 1: The evaluation of triple- index arrays of type Vijk = gijQik for three-body forces. 
The matrix with elements Qij is repeated along the third index k - forming the "pages" 
of three-index array. The same matrix is duplicated in the direction of the second index, 
filling the three-index array in another way. The total three-index array with elements 
Vijk is obtained in MATLAB by element-wise multiplication of the two "cubes" shown in 
the figure. The direction of the indeces is indicated. This figure illustrates the MATLAB 
treatment of eq. [251 



5.3 Evaluation of the forces from SW potential 

The total energy is sum over the triplets 





i>j>k 



but now the energies inside of the triplet must be specified. Since 



V{ri,rj,rk) = V{rj,ri,rk) 



V{rk,rj,ri) 
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With most reasonable interaction forms, V{ri, rj, r^) = Mjk+Mtij+^ki where the first index 
specifies the "vertex of the ordered triplet" , in the specific case discussed, the atom where 
the bond angle is measured. For such cases the summation can run over all individual 
cases and results in 

where now only factor 1/2 comes from the fact V^jk = Mkj? i-e. the vertex i is the same but 
the remaining pair (jk) can be counted twice. The evaluation of the forces is in principle 
identical, but can become quite different for cases where different degrees of simplification 
are possible. 

For The SW case, the interactions can be described as pair-separable. For a given 
Vijk one has the form 

% = 9irij)girii,)f{eij ■ e,fc) = g{rij)g{n\,) Q + e^^- ■ 6^^^ (26) 

Clearly, this form of interaction can be built as a three-index array, produced by element- 
wise multiplications of arrays constructed by duplication of matrices. Consider first a 
simpler form, where we keep just the funcions of distances, 

^ji = 9{rij)9{n\.) 

The total energy three-index array 

M = n Mr)) , (*, *, N)] ■ *7^ [{g{r)) , (*, iV, *)] 

The notation is schematic, and shows that the total array y^can be constructed as a 
product of two three index arrays each of which originates from replication of the original 
matrices. We introduce a symbol TZ [{A) ,(*,*, A^)] to denote the replication of matrix A 
along the dimension which contains the number, while the dimensions marked by asterix 
appear once. Construction of the two arrays is shown in figure [1] We have also used the 
operator A - *B for element-wise multiplication of two arrays of the same dimensions in 
analogy with the MATLAB notation. 

The importance of this is that also the gradient operations on such object will 
be obtained by operations on the original matrices, i.e. as Vmg{rij)- These are in form 
identical to the VmVij and so will be the Vrndik- These gradients thus need to be calculated 
only for the N x N matrix and replicated. Schematically we can write, following the above 
notation, 

[V] = n [(W^gir)) , (*, *, N)] ■ *7^ Mr)) , (*, A^, *)] 

+ n Mr)) ,(*,*, AT)] . *7^ [{Vmgir)) , (*, AT, *)] (27) 

Rewriting of these relations in the MATLAB notation will immediately produce the de- 
sired code. There will be a summation over two indices to produce the vector of forces 
from the three-index arrays. The evaluation of the gradients of the full SW-interaction 
follows in a straightforward manner the rules of chain derivative. The details and illustra- 
tion of the evaluation of the gradient of the cosine is given in the appendix. Implementing 
the result shown in the appendix, we can also there apply the discussed fact that the 
term consists of pair-separable terms. Thus all the operations can be also there per- 
formed on matrices of distances and duplicated to the necessary three-index objects in 
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precisely the same way as shown above for the distance functions. The only difference is 
that these are now Euclidian three dimensional vectors in addition to the particle-index 
dimensions. Thus the three-body forces can be constructed in a very similar way as the 
two-body forces. The computation is longer only due to the duplication of the matrices 
and element-wise products of the resulting matrices. The evaluation of the functional 
dependence does not change, since also here the computations are performed on matrices. 

5.4 Symmetry aspects of 3-body force computations 

The computations above can also be based on symmetry analysis. For a 3-body potential, 
each unique interacting triplet of particles has six different permutations contributing thus 
a total of six forces per index (particle). Since the Stillinger- Weber potential -formula \ijk, 
defined in eq. |23]- exhibits symmetry over the permutation of the second and third indices, 
i.e. Xijk = \ikj, so that forces on the k partilce satisfy 

^j(k) = Fi{k)j , F(^k)ij = F(k)ji , Fji(^k) = Fj(k)i (28) 

where the parenthesis denotes the particle coordinate of the gradient 

Fij^k) = -Vk{V,jk) (29) 

Relations in eq. [2H] demonstrate the symmetry of the problem and show that only three 
unique elements of the force array need to be computed for each particle in each triplet, as 
discussed also in other applications [6], [7]. Potential formula exhibits also the following 
antisymmetry due to the triplet's geometrical configuration; 

V,(l^,,fc) = -V,iV^,k) - Vfe(y,,,) (30) 

which can also be written as F{t)jk = — {^{j)k + Fij(k))- Thus only two unique forces per 
particle per unique triplet are to be evaluated. 

Computations of force arrays 

Evaluation of analytical derivatives (gradients) of the 3-body term of the given potential 
with respect to the third index (along the fourth dimension), gives the array, which is 
denoted F^, holding the following six unique (two on each particle) force elements per 
triplet: 



Fij{k) 


Fki{j) 


Fji{k) 


Fik{j) 


Fkj{i) 


Fjk(i) 



Those six unique force elements per triplet are what is needed to generate all the other 
elements by symmetry and antisymmetry exploitation. Transposing dimensions corre- 
sponding to the second and third indices of particles in F'^ array yields a new array F^, 
holding the other six force elements per triplet 



Fi{j)k 


Fk{i)j 


Fj(i)k 


Fi(k)j 


Fk(j)i 


Fj(k)i 



The third array F^ holding the last six force elements per triplet (cf. eq. [30]) 



F{i)jk = - 


'Fij^k) + 


= F{i)kj 


F(k)ij = - 




= F(k)ji 


F(j)ik = - 


[^m + Fj{i)k] 


= F^j)ki 
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is computed by the antisymmetry relation: 



F 



1 



(F^ + F^) 



Finally the total force on i-particle is a double sum over neighbors in triplets, i.e. 



This will be done in array form by addition of arrays F^, F^ [1 3 2 4] -permuted and F^ 
[1 4 3 2] -permuted (the permutation of array dimensions is done to place the particle i 
in the same position, see Appendix). The resulting 4-index array and then summing over 
the third and fourth index of results in 3 x array of the forces. 

Similar considerations and constructions can be done for the three-index arrays and 
their gradients in the case of Tersoff-type potentials in section O below, but there the 
symmetry and antisymmetry occur at two different stages of construction of the forces. 

5.5 Experiments with three-body forces 

In this section we have shown how one can use MATLAB functionality for evaluation of 
three body potentials and forces for three-body interactions. The operations on three- 
index arrays of the type discussed are "expensive" in computer time but rather straight- 
forward in terms of programming. The use of symmetries and more efficient build-up as 
discussed in sections 15.31 and 15.41 makes the computations much faster, but is much more 
error-prone. Here it can be useful to use a simple trick, formulate the computations first 
in a slow but easier programmed full-object version, even perhaps with the Kronnecker- 
deltas shown in eq. [15] implemented and explicitely summed, and then gradually introduce 
the more efficient computations in further versions. This offers possibilities for code de- 
velopment which are difficult to match in simplicity in other approaches. However, this 
approach offers even more benefits when investigating more complex interactions, as the 
Tersoff type discussed in the following section. 

6 TersofF-type many body potential 

A more complex many-body potentials quite extensively used are the Tersoff type poten- 
tials. They are constructed in a slightly different way, but the techniques discussed in 
previous section can be used with some modifications also here. The many body charac- 
ter is introduced via a pair interaction, but the parameters of the interaction depend on 
the positions of the neigbours and the same bond angles for triplets formed by the pair 
and a neighbour atom. We first follow the formulation usually found in literature and 
later translate this to our notation. Considering a Morse-like potential, the energy can 
be written as a sum of cite energies as follows: 





(31) 



where 




+ 




(32) 



13 



but where bij = b{rij; ri, r2, ra, ...fat), i.e. their value is dependent on all the other 
particles, in practice on those which are in the neighborhood, as assured by the cut-off 
function. The function fc{fij) is an aproperiate cutoff function, while fu and Ja are 
refered to as repulsive and attractive pair potentials respectively. 
The two potential functions of the pair-wise interactions are: 

fAin,) = -Se-^^^- ^^^^ 
The many-particle character in Tersoff's formula is introduced by the following: 

bij = (34) 

(l + /3"r^ij")27r 

Vij — Cijk 

Cijk = .Ur,k) gie^jk) e^i('---'--)' (35) 

'^^'^'^ = ^^i- i, ^^^^ 

" + [h- cos(%fc) ) 

where 9ijk is the angle between the and as defined above, and the set of parameters 
for carbon is given in [4J. In our work we used a Fermi-type cutoff function which has 
proven continuity and smoothness (infinitely differentiable) at any point, 

fc{r) 



r — fi ' 



e^—i + 1 

where i?j = i . 7, R2 = 2.0, = .5{R2 + Ri), = ^^(-Ra — i2i ), and it is defined and 
parametrized in correspondence to that introduced by Tersoff [5|, which has a different 
functional form but a very similar shape. Note that the same cut-off limits the pair 
interaction and the influence on the pair by the neighbors, the same fdi') appears both 
in eq. [32] and eq. [33 

To apply methods discussed in the previous section, we realize that the gradient 
operation can work on the various parts of the functional dependence. For simplicity, we 
omit the cut-off function in this discussion. The force on each particle m is calculated by 
applying the the gradient 

i j^i 

where we now re- write the above eqs. [31] to [3S] in a more transparent notation 

Vij= VR{rij) - VA{rij)B{Y,kg{rij)g{rik)d{eij ■ Gik)) (38) 



or 



where 



v., = V«(r,,) - VA{n,)B (Zj (39) 
Zij = ^ Qijk = E 9{rij)g{rik)d (e^^ ■ e^^) (40) 
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Applying the gradient, 

\Vm E E = ^ E E {^rnVRin,) - [ Va{t-.j)B (E, Qijk) ] ) (41) 

we recognize that the first term Vnlrij) is simply a two-body case discussed before. The 
really multi-body term is 

[ VAir,,)B (Zij) ] = V„ [ VAirij)B (J^k ^?(nj)(?(rik)rf (e., • e.^)) ] (42) 

The evaluation of the elements of matrices VA{rij) and Zij as well as the function B [Z^^ 
can be done in a straightforward way: construct the distance matrix as in eq [9l the unit 
vectors, and then construct the three-index array as in eq. in MATLAB notation 

Z=sum(Q , 3) ; 

produces the matrix Z. We can treat all of the terms as before, except the term 



VAin,)VmB ^E^(^ij)^?(^ik)c?(e., ■ 6,,^)^ = yA(r.,) V^S (Zy) = VAirij) ^^i^ 
but we recognize that there is only one new typ of object, 

E "^mQijk (43) 
k 

but the VmQijk has exactly the same form as the whole three-body contribution eq. [2^ in 
the SW case and thus can be evaluated and reduced to a matrix in the same way as in the 
SW case. The difference is that it must be combined with the other terms above. In spite 
of its complicated formal nature the Tersoff-type potentials still contain only quantities 
including triplets of particles and the construction of the code is only a slightly more 
complex modification of the method outlined in section 15.21 



7 Conclusion 

The discussed method, implementing multidimensional arrays along with their pre-defined 
operations, has proven to be a useful alternative to the standard loops aggregations in 
the vector-based MATLAB-like languages. It enables fast and powerful memory manip- 
ulations for interactive or exploratory work. In contrast to explicit loop programming 
it opens for speed-up and efficiency comparable to well optimized compiled arithmetic 
expressions in the framework of an interpreted code. Speed improvement is a main result 
of vectorized coding to which any numerical algorithm is very sensitive. Vectorization 
is known as performing operations on data stored as vectors, MATLAB stores matrices 
and arrays by columns in contiguous locations in the memory, and allows thus the use 
of vectorization on those objects. The elements of arrays can be selectively accessed in 
a vectorized fashion using so called logical indices. Mutidimensional logical arrays, as 
logical indices, have proven to be efficient and their use has resulted in a speed gain. 
Finally, the practicality, simplicity and generality of such a method to program empirical 
formulae for molecular dynamics purposes, make it effective to perform simulations on 
relatively small systems, up to 140 particles in a personal machine of 1 GB RAM, for 



15 



Tersoff's potential, a bigger system when using Stillinger-weber potential, and up to 1000 

particles if only a pair interactions are used. The described methodology can be applied 
both to investigation of aspects of molecular dynamics and to educational purposes. Fi- 
nally, programming in this method results in very short codes which are easily used and 
modified in interactive sessions. 



Appendices 

A Some aspects of MATLAB language 

Matrix(Array) manipulation in MATLAB 

Matrix(array) entering The lines starting with » are user input, the rest is 
matlab response. If input is ended by semicolon, no output follows. 



» 
A 

» 
B = 



A= [ 1 4 5 2.5 7 ] 



1.0000 
B= [ 1; 

1.0000 

4.0000 
2.5000 



4.0000 
4; 2.5] 



» 
C = 



1 

» D=[ 



C= [ 1 



D = 



1 0] 



2 ] 



1.0000 
3.0000 



4.0000 
7.0000 



5.0000 



3.0000 
6.6000 



2.5000 



7.0000 



5.0000 
2.0000 



2.0000 
2.0000 



The last example, matrix D shows multi-line possibility, the previous is a more 
compact matrix input. 

Accessing matrix(array) elements: Basic indexing 

Assume A is a random 6 by 6 matrix A=rand(6,6) 

A([l 2 3], 4) returns a column vector [ A(l,4) ; A(2,4) ; A(3,4) ] 

A(4, [1,1,1]) returns a row vector [ A(4,l) A(4,l) A(4,l) 1 

A([2,5] , [3,1]) returns a 2 by 2 matrix [ A(2,3) A(2,l) ; A(5,3) A(5,l) ] 

Matrix (Array) size, length and number of dimensions: 

The size of an array can be determined by the sizeO command as : 



[nrows ,ncols , npages, . . . .] = size(A) 



where 
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nrows=size(A, 1) , ncols=size(A,2) , npages=size(A,3) 

This function is of much importance when aiming elements-wise operations as it is 

the key for arrays' size matching. 

The ndims(A) function reports how many dimensions the array A has. Traihng 
singleton dimensions are ignored but singleton dimensions occurring before non- 
singleton dimensions are not. This means that 

a=5; c=a(l,l); B=zeros(5,5) ; B(5,5,l)=4; 

are valid expressions, even though they seemingly have non matching number of 
indices, but it is essential that the extra dimensions are "singletons", i.e. the extra 
dimensions have length one. 

An interesting construction: one can construct matrices and arrays expanding the 
singleton dimensions. 

a=5; B=a([l 1] , [1 1 1 1]) 
B = 

5 5 5 5 

5 5 5 5 

The length command gives the number of elements in the longest dimension, 
and is equivalent to max (size ()). It is also used with a row or a column vector. 

Matrix(Array) replication: 

One can use the MATLAB built-in script repmat, to replicate a matrix by the 
command: repmat (A , [m n] ) , which replicates the matrix A ( : , : ) m-times in the 
first dimension (rows), and n-times in the second dimension (columns). This is 
a block matrix construction, the block corresponding to the original matrix A is 
repeated in a m x n "supermatrix" . The operation repmat is not a built-in but it is 
a script which means a slow execution. An alternative efficient way to achieve the 
same replication is by : 

A( [1 : size(A, 1)] '*ones(l,m) , [1 : size(A,2)] '*ones(l,n)) which is an applica- 
tion of basic indexing techniques. 

In this work we can use another aspect of repmat () script, for building multidi- 
mensional arrays from matrices. Given a matrix A, its elements can be accessed by 
A( : , : ) , but that is equivalent to A( : , : , 1) , but also A( : , : , 1 , 1) , this is referred 
to as "trailing singleton dimensions". Using this, one can replicate the traihng di- 
mensions. From an NxN matrix we can construct a three- index NxNxN array 
by B=repmat(A, [1,1,N]); In the text we discuss a more efficient equivalent 
B=A(: , : ,ones(l,N) ) 

Transpose and its generalization permute: 

For a matrix (2D array) non-conjugate transpose is given by transpose (A) and 

it is equivalent to a standard matrix operation A . ' . The normal Hcrmitc type 
conjugation is denoted by A' and has an equivalent function ctranspose(A). The 
A. '-notation is somewhat surprising, since the dot usually denotes element- wise 
counterpart of the undotted operator (as e.g. A. "2, or A .* B etc). 
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The generalizations of transpose-function to multi-dimensional arrays are given as 
permute (A, ORDER) and ipermute(A,ORDER) where ORDER means an array 
of reordered indeces. When ORDER is [1 2 3 4] for 4-index array, it means identity. 

transpose(A) and permute(A, [2 1]) give the same results as (A.'). The two 
permute-functions are related as shown in the example 

H=permute (G , [3 , 1 , 2] ) , G=ipermute (H , [3 , 1 , 2] ) 

If size (G) is e.g. 3 4 2 then size (H) is 2 3 4 

Clearly, the inverse permute is introduced for convenience, the same result could be 
achieved by using the permute with appropriate order of dimensions. 

Logical array addressing and predefined functions: 

For a given matrix or array G the function logical (G) returns an equal size array 
where there are logical ones (true) in elements which were nonzero and logical zeros 
(false) where the original array had zeros. Thus 

Lg=logical(G) and Lg=G>0 both return the same matrix of logical ones and zeros. 

The extreme usefulness of these functions is that when we write 

G(Lg)= 1 ./ G(Lg) the calculation will be performed only for the logically true 
addresses, thus no division by zero will happen. Or, we can "shrink" a range of 
values; for example all values below a threshold can be given the threshold value 
and all values above a chosen top value can be given the top value. As an example, 
we take a matrix which should be modified by replacing all values smaller than 2 
by 2 and all larger than 5 by 5. This can be done by the following commands 

B=A; B(A>5) = A(A>5)*0+5; B(A<2) = A(A<2)*0+2 

The original matrix and the resulting one are 

A=[1536 B= 2535 

7311 5322 
9341 5342 
3421 3422 
9231] 5232 

There are also available functions analogous to ones( ) and zeros ( ). As 
an example, false (5,2,2) returns the same logical array as zeros (5 , 2 , 2) >0 and 
true (5,2,2) is the same as e.g. logical (ones (5,2,2)) 

The use of logical addressing can shorten some large array manipulations and 
calculations, but also consumes some processing time. 

One can test the gains and losses by filling a large array and doing the element- 
wise operations on the whole array and the logically addressed part of the array. 

In this example we show the array of ten million numbers, where only 100 are 
nonzero. 

N = 10000000; 

octave > tic; C=zeros(N, 1) ; C(100 : 200)=1 ; D=exp(C); toe 
Elapsed time is 0.601 seconds. 

octave > tic;C=zeros(N,l) ;C(100:200)=1; D(C>0)=exp(C(C>0) ) ;toc 
Elapsed time is 0.394 seconds. 
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In the second line we have calculated only the nonzero cases using the discussed 
logical addressing. We observe that there is a certain small gain in time, but far 
under the factor 10^ expected from the number of calculations. It should also be 
mentioned that the return times are not reproduceable, they depend on the over all 
state of the computer memory. 

Speed improvement in Matlab 

Array Preallocation : 

Matlabs matrix variables dynamically augment rows and columns and the matrix is 
automatically resized. Internally, the matrix data memory must be reallocated with 
larger size. If a matrix is resized repeatedly like within a for loop, this overhead 
becomes noticeable. When a loop is of a must demand, as the inevitable time-loop 
in MD programs, frequent reallocations are avoided by preallocating the array with 
one of the zeros (), onesO, randO, trueO and false () commands. 

Vectorization of MATLAB 's computation and logical addressing: Vectorized 
computations which mean replacing parallel operations, i.e., simultaneous execution 
of those operations in loops, with vector operations, improves the speed often ten- 
fold. For that purpose, most standard functions in Matlab are vectorized, i.e., they 
operate on an array as if the function was applied individually to every element. 
Logical addressing operations are also vectorized, the use of logical addressing may 
lead to vectorized computations only for elements where the indices are true, thus 
reducing execution time, in some cases substantially. 

B Gradient of the cosine of the bond angle for three 
body forces 

Let a, b and c mark three atoms, the b-atom is at the vertex of the angle 6abc', 
Define the following quantities: 

Tfec = — Tb = fbc^bc] ^ba = Ta " ^fe = ^ba^ba] COS Babe = ^ba ' ^bc] 



a 




c 



e 



be 



b 



be 



Gradient with respect to the end 



c 




(44) 



Gradient with respect to the end a 



Va cos 9, 



1 



[etc - efe„ cos 6 abc] 



(45) 



abc 



gradient with respect to point 'b' : 

Vft cos 6 abc = — Gba- 

Tbc 



etc COS 9abc - e;,^ cos Babe 

Tba 



(46) 
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The gradients satisfy the following relation 



Vh cos 9, 



abc 



- Va COS dabc - Vc COS d, 



abc 



(47) 



which simply confirms the fact that moving the vertex point by a vector Ar is the same 
as moving each of the end points a and c simultaneously by the vector — Ar. 
Derivation: 

In view of eq. HT] it is necessary to evaluate just one endpoint, the second follows 
from symmetry and the vertex is given by relation eq. HT] 
The calculation is done by using 



and 



Vj-(r ■ const) ^ const 



1 _ _^ _ _^ 



Applying these two to the chain derivatives in the gradient as 



^bc ^ba 

nc na 



nc na 





1 


^ba 


^bc 


( ^bc 


rba 


rL na 


nc 


na 


nc 


\nc 


na 



(48) 
(49) 
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